/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "distributedTriSurfaceMesh.H"
#include "mapDistribute.H"
#include "Random.H"
#include "addToRunTimeSelectionTable.H"
#include "triangleFuncs.H"
#include "matchPoints.H"
#include "globalIndex.H"
#include "foamTime.H"

#include "IFstream.H"
#include "decompositionMethod.H"
#include "vectorList.H"
#include "PackedBoolList.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
	defineTypeNameAndDebug(distributedTriSurfaceMesh, 0);
	addToRunTimeSelectionTable
	(
		searchableSurface,
		distributedTriSurfaceMesh,
		dict
	);
}


template<>
const char*
Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>::names[] =
{
	"follow",
	"independent",
	"frozen"
};

const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
	Foam::distributedTriSurfaceMesh::distributionTypeNames_;


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

// Read my additional data from the dictionary
bool Foam::distributedTriSurfaceMesh::read()
{
	// Get bb of all domains.
	procBb_.setSize(Pstream::nProcs());

	procBb_[Pstream::myProcNo()] = List<treeBoundBox>(dict_.lookup("bounds"));
	Pstream::gatherList(procBb_);
	Pstream::scatterList(procBb_);

	// Distribution type
	distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));

	// Merge distance
	mergeDist_ = readScalar(dict_.lookup("mergeDistance"));

	return true;
}


// Is segment fully local?
bool Foam::distributedTriSurfaceMesh::isLocal
(
	const List<treeBoundBox>& myBbs,
	const point& start,
	const point& end
)
{
	forAll(myBbs, bbI)
	{
		if (myBbs[bbI].contains(start) && myBbs[bbI].contains(end))
		{
			return true;
		}
	}
	return false;
}


//void Foam::distributedTriSurfaceMesh::splitSegment
//(
//    const label segmentI,
//    const point& start,
//    const point& end,
//    const treeBoundBox& bb,
//
//    DynamicList<segment>& allSegments,
//    dynamicLabelList& allSegmentMap,
//    dynamicLabelList sendMap
//) const
//{
//    // Work points
//    point clipPt0, clipPt1;
//
//    if (bb.contains(start))
//    {
//        // start within, trim end to bb
//        bool clipped = bb.intersects(end, start, clipPt0);
//
//        if (clipped)
//        {
//            // segment from start to clippedStart passes
//            // through proc.
//            sendMap[procI].append(allSegments.size());
//            allSegmentMap.append(segmentI);
//            allSegments.append(segment(start, clipPt0));
//        }
//    }
//    else if (bb.contains(end))
//    {
//        // end within, trim start to bb
//        bool clipped = bb.intersects(start, end, clipPt0);
//
//        if (clipped)
//        {
//            sendMap[procI].append(allSegments.size());
//            allSegmentMap.append(segmentI);
//            allSegments.append(segment(clipPt0, end));
//        }
//    }
//    else
//    {
//        // trim both
//        bool clippedStart = bb.intersects(start, end, clipPt0);
//
//        if (clippedStart)
//        {
//            bool clippedEnd = bb.intersects(end, clipPt0, clipPt1);
//
//            if (clippedEnd)
//            {
//                // middle part of segment passes through proc.
//                sendMap[procI].append(allSegments.size());
//                allSegmentMap.append(segmentI);
//                allSegments.append(segment(clipPt0, clipPt1));
//            }
//        }
//    }
//}


void Foam::distributedTriSurfaceMesh::distributeSegment
(
	const label segmentI,
	const point& start,
	const point& end,

	DynamicList<segment>& allSegments,
	dynamicLabelList& allSegmentMap,
	List<dynamicLabelList >& sendMap
) const
{
	// Work points
	point clipPt;


	// 1. Fully local already handled outside. Note: retest is cheap.
	if (isLocal(procBb_[Pstream::myProcNo()], start, end))
	{
		return;
	}


	// 2. If fully inside one other processor, then only need to send
	// to that one processor even if it intersects another. Rare occurrence
	// but cheap to test.
	forAll(procBb_, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			const List<treeBoundBox>& bbs = procBb_[procI];

			if (isLocal(bbs, start, end))
			{
				sendMap[procI].append(allSegments.size());
				allSegmentMap.append(segmentI);
				allSegments.append(segment(start, end));
				return;
			}
		}
	}


	// 3. If not contained in single processor send to all intersecting
	// processors.
	forAll(procBb_, procI)
	{
		const List<treeBoundBox>& bbs = procBb_[procI];

		forAll(bbs, bbI)
		{
			const treeBoundBox& bb = bbs[bbI];

			// Scheme a: any processor that intersects the segment gets
			// the segment.

			if (bb.intersects(start, end, clipPt))
			{
				sendMap[procI].append(allSegments.size());
				allSegmentMap.append(segmentI);
				allSegments.append(segment(start, end));
			}

			// Alternative: any processor only gets clipped bit of
			// segment. This gives small problems with additional
			// truncation errors.
			//splitSegment
			//(
			//    segmentI,
			//    start,
			//    end,
			//    bb,
			//
			//    allSegments,
			//    allSegmentMap,
			//   sendMap[procI]
			//);
		}
	}
}


Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::distributeSegments
(
	const pointField& start,
	const pointField& end,

	List<segment>& allSegments,
	labelList& allSegmentMap
) const
{
	// Determine send map
	// ~~~~~~~~~~~~~~~~~~

	labelListList sendMap(Pstream::nProcs());

	{
		// Since intersection test is quite expensive compared to memory
		// allocation we use DynamicList to immediately store the segment
		// in the correct bin.

		// Segments to test
		DynamicList<segment> dynAllSegments(start.size());
		// Original index of segment
		dynamicLabelList dynAllSegmentMap(start.size());
		// Per processor indices into allSegments to send
		List<dynamicLabelList > dynSendMap(Pstream::nProcs());

		forAll(start, segmentI)
		{
			distributeSegment
			(
				segmentI,
				start[segmentI],
				end[segmentI],

				dynAllSegments,
				dynAllSegmentMap,
				dynSendMap
			);
		}

		// Convert dynamicList to labelList
		sendMap.setSize(Pstream::nProcs());
		forAll(sendMap, procI)
		{
			dynSendMap[procI].shrink();
			sendMap[procI].transfer(dynSendMap[procI]);
		}

		allSegments.transfer(dynAllSegments.shrink());
		allSegmentMap.transfer(dynAllSegmentMap.shrink());
	}


	// Send over how many I need to receive.
	labelListList sendSizes(Pstream::nProcs());
	sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
	forAll(sendMap, procI)
	{
		sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
	}
	Pstream::gatherList(sendSizes);
	Pstream::scatterList(sendSizes);


	// Determine order of receiving
	labelListList constructMap(Pstream::nProcs());

	// My local segments first
	constructMap[Pstream::myProcNo()] = identity
	(
		sendMap[Pstream::myProcNo()].size()
	);

	label segmentI = constructMap[Pstream::myProcNo()].size();
	forAll(constructMap, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			// What I need to receive is what other processor is sending to me.
			label nRecv = sendSizes[procI][Pstream::myProcNo()];
			constructMap[procI].setSize(nRecv);

			for (label i = 0; i < nRecv; i++)
			{
				constructMap[procI][i] = segmentI++;
			}
		}
	}

	return autoPtr<mapDistribute>
	(
		new mapDistribute
		(
			segmentI,       // size after construction
			sendMap,
			constructMap,
			true            // reuse storage
		)
	);
}


void Foam::distributedTriSurfaceMesh::findLine
(
	const bool nearestIntersection,
	const pointField& start,
	const pointField& end,
	List<pointIndexHit>& info
) const
{
	const indexedOctree<treeDataTriSurface>& octree = tree();

	// Important:force synchronised construction of indexing
	const globalIndex& triIndexer = globalTris();

	// Initialise
	info.setSize(start.size());
	forAll(info, i)
	{
		info[i].setMiss();
	}


	// Do any local queries
	// ~~~~~~~~~~~~~~~~~~~~

	label nLocal = 0;

	forAll(start, i)
	{
		if (isLocal(procBb_[Pstream::myProcNo()], start[i], end[i]))
		{
			if (nearestIntersection)
			{
				info[i] = octree.findLine(start[i], end[i]);
			}
			else
			{
				info[i] = octree.findLineAny(start[i], end[i]);
			}

			if (info[i].hit())
			{
				info[i].setIndex(triIndexer.toGlobal(info[i].index()));
			}
			nLocal++;
		}
	}


	if
	(
		Pstream::parRun()
	 && (
			returnReduce(nLocal, sumOp<label>())
		  < returnReduce(start.size(), sumOp<label>())
		)
	)
	{
		// Not all can be resolved locally. Build segments and map, send over
		// segments, do intersections, send back and merge.


		// Construct queries (segments)
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		// Segments to test
		List<segment> allSegments(start.size());
		// Original index of segment
		labelList allSegmentMap(start.size());

		const autoPtr<mapDistribute> mapPtr
		(
			distributeSegments
			(
				start,
				end,
				allSegments,
				allSegmentMap
			)
		);
		const mapDistribute& map = mapPtr();

		label nOldAllSegments = allSegments.size();


		// Exchange the segments
		// ~~~~~~~~~~~~~~~~~~~~~

		map.distribute
		(
			Pstream::nonBlocking,   //Pstream::scheduled,
			List<labelPair>(0),     //map.schedule(),
			map.constructSize(),
			map.subMap(),           // what to send
			map.constructMap(),     // what to receive
			allSegments
		);


		// Do tests I need to do
		// ~~~~~~~~~~~~~~~~~~~~~

		// Intersections
		List<pointIndexHit> intersections(allSegments.size());

		forAll(allSegments, i)
		{
			if (nearestIntersection)
			{
				intersections[i] = octree.findLine
				(
					allSegments[i].first(),
					allSegments[i].second()
				);
			}
			else
			{
				intersections[i] = octree.findLineAny
				(
					allSegments[i].first(),
					allSegments[i].second()
				);
			}

			// Convert triangle index to global numbering
			if (intersections[i].hit())
			{
				intersections[i].setIndex
				(
					triIndexer.toGlobal(intersections[i].index())
				);
			}
		}


		// Exchange the intersections (opposite to segments)
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		map.distribute
		(
			//Pstream::scheduled,
			//map.schedule            // Note reverse schedule
			//(
			//    map.constructMap(),
			//    map.subMap()
			//),
			Pstream::nonBlocking,
			List<labelPair>(0),
			nOldAllSegments,
			map.constructMap(),     // what to send
			map.subMap(),           // what to receive
			intersections
		);


		// Extract the hits
		// ~~~~~~~~~~~~~~~~

		forAll(intersections, i)
		{
			const pointIndexHit& allInfo = intersections[i];
			label segmentI = allSegmentMap[i];
			pointIndexHit& hitInfo = info[segmentI];

			if (allInfo.hit())
			{
				if (!hitInfo.hit())
				{
					// No intersection yet so take this one
					hitInfo = allInfo;
				}
				else if (nearestIntersection)
				{
					// Nearest intersection
					if
					(
					    magSqr(allInfo.hitPoint()-start[segmentI])
					  < magSqr(hitInfo.hitPoint()-start[segmentI])
					)
					{
					    hitInfo = allInfo;
					}
				}
			}
		}
	}
}


// Exchanges indices to the processor they come from.
// - calculates exchange map
// - uses map to calculate local triangle index
Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
	const List<pointIndexHit>& info,
	labelList& triangleIndex
) const
{
	triangleIndex.setSize(info.size());

	const globalIndex& triIndexer = globalTris();


	// Determine send map
	// ~~~~~~~~~~~~~~~~~~

	// Since determining which processor the query should go to is
	// cheap we do a multi-pass algorithm to save some memory temporarily.

	// 1. Count
	labelList nSend(Pstream::nProcs(), 0);

	forAll(info, i)
	{
		if (info[i].hit())
		{
			label procI = triIndexer.whichProcID(info[i].index());
			nSend[procI]++;
		}
	}

	// 2. Size sendMap
	labelListList sendMap(Pstream::nProcs());
	forAll(nSend, procI)
	{
		sendMap[procI].setSize(nSend[procI]);
		nSend[procI] = 0;
	}

	// 3. Fill sendMap
	forAll(info, i)
	{
		if (info[i].hit())
		{
			label procI = triIndexer.whichProcID(info[i].index());
			triangleIndex[i] = triIndexer.toLocal(procI, info[i].index());
			sendMap[procI][nSend[procI]++] = i;
		}
		else
		{
			triangleIndex[i] = -1;
		}
	}


	// Send over how many I need to receive
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	labelListList sendSizes(Pstream::nProcs());
	sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
	forAll(sendMap, procI)
	{
		sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
	}
	Pstream::gatherList(sendSizes);
	Pstream::scatterList(sendSizes);


	// Determine receive map
	// ~~~~~~~~~~~~~~~~~~~~~

	labelListList constructMap(Pstream::nProcs());

	// My local segments first
	constructMap[Pstream::myProcNo()] = identity
	(
		sendMap[Pstream::myProcNo()].size()
	);

	label segmentI = constructMap[Pstream::myProcNo()].size();
	forAll(constructMap, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			// What I need to receive is what other processor is sending to me.
			label nRecv = sendSizes[procI][Pstream::myProcNo()];
			constructMap[procI].setSize(nRecv);

			for (label i = 0; i < nRecv; i++)
			{
				constructMap[procI][i] = segmentI++;
			}
		}
	}


	// Pack into distribution map
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~

	autoPtr<mapDistribute> mapPtr
	(
		new mapDistribute
		(
			segmentI,       // size after construction
			sendMap,
			constructMap,
			true            // reuse storage
		)
	);
	const mapDistribute& map = mapPtr();


	// Send over queries
	// ~~~~~~~~~~~~~~~~~

	map.distribute
	(
		//Pstream::scheduled,
		//map.schedule(),
		Pstream::nonBlocking,
		List<labelPair>(0),
		map.constructSize(),
		map.subMap(),           // what to send
		map.constructMap(),     // what to receive
		triangleIndex
	);


	return mapPtr;
}


Foam::label Foam::distributedTriSurfaceMesh::calcOverlappingProcs
(
	const point& centre,
	const scalar radiusSqr,
	boolList& overlaps
) const
{
	overlaps = false;
	label nOverlaps = 0;

	forAll(procBb_, procI)
	{
		const List<treeBoundBox>& bbs = procBb_[procI];

		forAll(bbs, bbI)
		{
			if (bbs[bbI].overlaps(centre, radiusSqr))
			{
				overlaps[procI] = true;
				nOverlaps++;
				break;
			}
		}
	}
	return nOverlaps;
}


// Generate queries for parallel distance calculation
// - calculates exchange map
// - uses map to exchange points and radius
Foam::autoPtr<Foam::mapDistribute>
Foam::distributedTriSurfaceMesh::calcLocalQueries
(
	const pointField& centres,
	const scalarField& radiusSqr,

	pointField& allCentres,
	scalarField& allRadiusSqr,
	labelList& allSegmentMap
) const
{
	// Determine queries
	// ~~~~~~~~~~~~~~~~~

	labelListList sendMap(Pstream::nProcs());

	{
		// Queries
		DynamicList<point> dynAllCentres(centres.size());
		DynamicList<scalar> dynAllRadiusSqr(centres.size());
		// Original index of segment
		dynamicLabelList dynAllSegmentMap(centres.size());
		// Per processor indices into allSegments to send
		List<dynamicLabelList > dynSendMap(Pstream::nProcs());

		// Work array - whether processor bb overlaps the bounding sphere.
		boolList procBbOverlaps(Pstream::nProcs());

		forAll(centres, centreI)
		{
			// Find the processor this sample+radius overlaps.
			calcOverlappingProcs
			(
				centres[centreI],
				radiusSqr[centreI],
				procBbOverlaps
			);

			forAll(procBbOverlaps, procI)
			{
				if (procI != Pstream::myProcNo() && procBbOverlaps[procI])
				{
					dynSendMap[procI].append(dynAllCentres.size());
					dynAllSegmentMap.append(centreI);
					dynAllCentres.append(centres[centreI]);
					dynAllRadiusSqr.append(radiusSqr[centreI]);
				}
			}
		}

		// Convert dynamicList to labelList
		sendMap.setSize(Pstream::nProcs());
		forAll(sendMap, procI)
		{
			dynSendMap[procI].shrink();
			sendMap[procI].transfer(dynSendMap[procI]);
		}

		allCentres.transfer(dynAllCentres.shrink());
		allRadiusSqr.transfer(dynAllRadiusSqr.shrink());
		allSegmentMap.transfer(dynAllSegmentMap.shrink());
	}


	// Send over how many I need to receive.
	labelListList sendSizes(Pstream::nProcs());
	sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
	forAll(sendMap, procI)
	{
		sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
	}
	Pstream::gatherList(sendSizes);
	Pstream::scatterList(sendSizes);


	// Determine order of receiving
	labelListList constructMap(Pstream::nProcs());

	// My local segments first
	constructMap[Pstream::myProcNo()] = identity
	(
		sendMap[Pstream::myProcNo()].size()
	);

	label segmentI = constructMap[Pstream::myProcNo()].size();
	forAll(constructMap, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			// What I need to receive is what other processor is sending to me.
			label nRecv = sendSizes[procI][Pstream::myProcNo()];
			constructMap[procI].setSize(nRecv);

			for (label i = 0; i < nRecv; i++)
			{
				constructMap[procI][i] = segmentI++;
			}
		}
	}

	autoPtr<mapDistribute> mapPtr
	(
		new mapDistribute
		(
			segmentI,       // size after construction
			sendMap,
			constructMap,
			true            // reuse storage
		)
	);
	return mapPtr;
}


// Find bounding boxes that guarantee a more or less uniform distribution
// of triangles. Decomposition in here is only used to get the bounding
// boxes, actual decomposition is done later on.
// Returns a per processor a list of bounding boxes that most accurately
// describe the shape. For now just a single bounding box per processor but
// optimisation might be to determine a better fitting shape.
Foam::List<Foam::List<Foam::treeBoundBox> >
Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
(
	const triSurface& s
)
{
	if (!decomposer_.valid())
	{
		// Use current decomposer.
		// Note: or always use hierarchical?
		IOdictionary decomposeDict
		(
			IOobject
			(
				"decomposeParDict",
				searchableSurface::time().system(),
				searchableSurface::time(),
				IOobject::MUST_READ_IF_MODIFIED,
				IOobject::NO_WRITE,
				false
			)
		);
		decomposer_ = decompositionMethod::New(decomposeDict);

		if (!decomposer_().parallelAware())
		{
			FatalErrorIn
			(
				"distributedTriSurfaceMesh::independentlyDistributedBbs"
				"(const triSurface&)"
			)   << "The decomposition method " << decomposer_().typeName
				<< " does not decompose in parallel."
				<< " Please choose one that does." << exit(FatalError);
		}
	}

	// Do decomposition according to triangle centre
	pointField triCentres(s.size());
	forAll (s, triI)
	{
		triCentres[triI] = s[triI].centre(s.points());
	}

	// Do the actual decomposition
	labelList distribution(decomposer_->decompose(triCentres));

	// Find bounding box for all triangles on new distribution.

	// Initialise to inverted box (VGREAT, -VGREAT)
	List<List<treeBoundBox> > bbs(Pstream::nProcs());
	forAll(bbs, procI)
	{
		bbs[procI].setSize(1);
		//bbs[procI][0] = boundBox::invertedBox;
		bbs[procI][0].min() = point( VGREAT,  VGREAT,  VGREAT);
		bbs[procI][0].max() = point(-VGREAT, -VGREAT, -VGREAT);
	}

	forAll (s, triI)
	{
		point& bbMin = bbs[distribution[triI]][0].min();
		point& bbMax = bbs[distribution[triI]][0].max();

		const labelledTri& f = s[triI];
		const point& p0 = s.points()[f[0]];
		const point& p1 = s.points()[f[1]];
		const point& p2 = s.points()[f[2]];

		bbMin = min(bbMin, p0);
		bbMin = min(bbMin, p1);
		bbMin = min(bbMin, p2);

		bbMax = max(bbMax, p0);
		bbMax = max(bbMax, p1);
		bbMax = max(bbMax, p2);
	}

	// Now combine for all processors and convert to correct format.
	forAll(bbs, procI)
	{
		forAll(bbs[procI], i)
		{
			reduce(bbs[procI][i].min(), minOp<point>());
			reduce(bbs[procI][i].max(), maxOp<point>());
		}
	}
	return bbs;
}


// Does any part of triangle overlap bb.
bool Foam::distributedTriSurfaceMesh::overlaps
(
	const List<treeBoundBox>& bbs,
	const point& p0,
	const point& p1,
	const point& p2
)
{
	forAll(bbs, bbI)
	{
		const treeBoundBox& bb = bbs[bbI];

		treeBoundBox triBb(p0, p0);
		triBb.min() = min(triBb.min(), p1);
		triBb.min() = min(triBb.min(), p2);

		triBb.max() = max(triBb.max(), p1);
		triBb.max() = max(triBb.max(), p2);

		//- Exact test of triangle intersecting bb

		// Quick rejection. If whole bounding box of tri is outside cubeBb then
		// there will be no intersection.
		if (bb.overlaps(triBb))
		{
			// Check if one or more triangle point inside
			if (bb.contains(p0) || bb.contains(p1) || bb.contains(p2))
			{
				// One or more points inside
				return true;
			}

			// Now we have the difficult case: all points are outside but
			// connecting edges might go through cube. Use fast intersection
			// of bounding box.
			bool intersect = triangleFuncs::intersectBb(p0, p1, p2, bb);

			if (intersect)
			{
				return true;
			}
		}
	}
	return false;
}


void Foam::distributedTriSurfaceMesh::subsetMeshMap
(
	const triSurface& s,
	const boolList& include,
	const label nIncluded,
	labelList& newToOldPoints,
	labelList& oldToNewPoints,
	labelList& newToOldFaces
)
{
	newToOldFaces.setSize(nIncluded);
	newToOldPoints.setSize(s.points().size());
	oldToNewPoints.setSize(s.points().size());
	oldToNewPoints = -1;
	{
		label faceI = 0;
		label pointI = 0;

		forAll(include, oldFacei)
		{
			if (include[oldFacei])
			{
				// Store new faces compact
				newToOldFaces[faceI++] = oldFacei;

				// Renumber labels for triangle
				const labelledTri& tri = s[oldFacei];

				forAll(tri, fp)
				{
					label oldPointI = tri[fp];

					if (oldToNewPoints[oldPointI] == -1)
					{
					    oldToNewPoints[oldPointI] = pointI;
					    newToOldPoints[pointI++] = oldPointI;
					}
				}
			}
		}
		newToOldPoints.setSize(pointI);
	}
}


Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
	const triSurface& s,
	const labelList& newToOldPoints,
	const labelList& oldToNewPoints,
	const labelList& newToOldFaces
)
{
	// Extract points
	pointField newPoints(newToOldPoints.size());
	forAll(newToOldPoints, i)
	{
		newPoints[i] = s.points()[newToOldPoints[i]];
	}
	// Extract faces
	List<labelledTri> newTriangles(newToOldFaces.size());

	forAll(newToOldFaces, i)
	{
		// Get old vertex labels
		const labelledTri& tri = s[newToOldFaces[i]];

		newTriangles[i][0] = oldToNewPoints[tri[0]];
		newTriangles[i][1] = oldToNewPoints[tri[1]];
		newTriangles[i][2] = oldToNewPoints[tri[2]];
		newTriangles[i].region() = tri.region();
	}

	// Reuse storage.
	return triSurface(newTriangles, s.patches(), newPoints, true);
}


Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
	const triSurface& s,
	const boolList& include,
	labelList& newToOldPoints,
	labelList& newToOldFaces
)
{
	label n = 0;

	forAll(include, i)
	{
		if (include[i])
		{
			n++;
		}
	}

	labelList oldToNewPoints;
	subsetMeshMap
	(
		s,
		include,
		n,
		newToOldPoints,
		oldToNewPoints,
		newToOldFaces
	);

	return subsetMesh
	(
		s,
		newToOldPoints,
		oldToNewPoints,
		newToOldFaces
	);
}


Foam::triSurface Foam::distributedTriSurfaceMesh::subsetMesh
(
	const triSurface& s,
	const labelList& newToOldFaces,
	labelList& newToOldPoints
)
{
	const boolList include
	(
		createWithValues<boolList>
		(
			s.size(),
			false,
			newToOldFaces,
			true
		)
	);

	newToOldPoints.setSize(s.points().size());
	labelList oldToNewPoints(s.points().size(), -1);
	{
		label pointI = 0;

		forAll(include, oldFacei)
		{
			if (include[oldFacei])
			{
				// Renumber labels for triangle
				const labelledTri& tri = s[oldFacei];

				forAll(tri, fp)
				{
					label oldPointI = tri[fp];

					if (oldToNewPoints[oldPointI] == -1)
					{
					    oldToNewPoints[oldPointI] = pointI;
					    newToOldPoints[pointI++] = oldPointI;
					}
				}
			}
		}
		newToOldPoints.setSize(pointI);
	}

	return subsetMesh
	(
		s,
		newToOldPoints,
		oldToNewPoints,
		newToOldFaces
	);
}


Foam::label Foam::distributedTriSurfaceMesh::findTriangle
(
	const List<labelledTri>& allFaces,
	const labelListList& allPointFaces,
	const labelledTri& otherF
)
{
	// allFaces connected to otherF[0]
	const labelList& pFaces = allPointFaces[otherF[0]];

	forAll(pFaces, i)
	{
		const labelledTri& f = allFaces[pFaces[i]];

		if (f.region() == otherF.region())
		{
			// Find index of otherF[0]
			label fp0 = findIndex(f, otherF[0]);
			// Check rest of triangle in same order
			label fp1 = f.fcIndex(fp0);
			label fp2 = f.fcIndex(fp1);

			if (f[fp1] == otherF[1] && f[fp2] == otherF[2])
			{
				return pFaces[i];
			}
		}
	}
	return -1;
}


// Merge into allSurf.
void Foam::distributedTriSurfaceMesh::merge
(
	const scalar mergeDist,
	const List<labelledTri>& subTris,
	const pointField& subPoints,

	List<labelledTri>& allTris,
	pointField& allPoints,

	labelList& faceConstructMap,
	labelList& pointConstructMap
)
{
	labelList subToAll;
	matchPoints
	(
		subPoints,
		allPoints,
		scalarField(subPoints.size(), mergeDist),   // match distance
		false,                                      // verbose
		pointConstructMap
	);

	label nOldAllPoints = allPoints.size();


	// Add all unmatched points
	// ~~~~~~~~~~~~~~~~~~~~~~~~

	label allPointI = nOldAllPoints;
	forAll(pointConstructMap, pointI)
	{
		if (pointConstructMap[pointI] == -1)
		{
			pointConstructMap[pointI] = allPointI++;
		}
	}

	if (allPointI > nOldAllPoints)
	{
		allPoints.setSize(allPointI);

		forAll(pointConstructMap, pointI)
		{
			if (pointConstructMap[pointI] >= nOldAllPoints)
			{
				allPoints[pointConstructMap[pointI]] = subPoints[pointI];
			}
		}
	}


	// To check whether triangles are same we use pointFaces.
	labelListList allPointFaces;
	invertManyToMany(nOldAllPoints, allTris, allPointFaces);


	// Add all unmatched triangles
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~

	label allTriI = allTris.size();
	allTris.setSize(allTriI + subTris.size());

	faceConstructMap.setSize(subTris.size());

	forAll(subTris, triI)
	{
		const labelledTri& subTri = subTris[triI];

		// Get triangle in new numbering
		labelledTri mappedTri
		(
			pointConstructMap[subTri[0]],
			pointConstructMap[subTri[1]],
			pointConstructMap[subTri[2]],
			subTri.region()
		);


		// Check if all points were already in surface
		bool fullMatch = true;

		forAll(mappedTri, fp)
		{
			if (mappedTri[fp] >= nOldAllPoints)
			{
				fullMatch = false;
				break;
			}
		}

		if (fullMatch)
		{
			// All three points are mapped to old points. See if same
			// triangle.
			label i = findTriangle
			(
				allTris,
				allPointFaces,
				mappedTri
			);

			if (i == -1)
			{
				// Add
				faceConstructMap[triI] = allTriI;
				allTris[allTriI] = mappedTri;
				allTriI++;
			}
			else
			{
				faceConstructMap[triI] = i;
			}
		}
		else
		{
			// Add
			faceConstructMap[triI] = allTriI;
			allTris[allTriI] = mappedTri;
			allTriI++;
		}
	}
	allTris.setSize(allTriI);
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
(
	const IOobject& io,
	const triSurface& s,
	const dictionary& dict
)
:
	triSurfaceMesh(io, s),
	dict_
	(
		IOobject
		(
			searchableSurface::name() + "Dict",
			searchableSurface::instance(),
			searchableSurface::local(),
			searchableSurface::db(),
			searchableSurface::NO_READ,
			searchableSurface::writeOpt(),
			searchableSurface::registerObject()
		),
		dict
	)
{
	read();

	if (debug)
	{
		Info<< "Constructed from triSurface:" << endl;
		writeStats(Info);

		labelList nTris(Pstream::nProcs());
		nTris[Pstream::myProcNo()] = triSurface::size();
		Pstream::gatherList(nTris);
		Pstream::scatterList(nTris);

		Info<< endl<< "\tproc\ttris\tbb" << endl;
		forAll(nTris, procI)
		{
			Info<< '\t' << procI << '\t' << nTris[procI]
				 << '\t' << procBb_[procI] << endl;
		}
		Info<< endl;
	}
}


Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
:
	//triSurfaceMesh(io),
	triSurfaceMesh
	(
		IOobject
		(
			io.name(),
			io.time().findInstance(io.local(), word::null),
			io.local(),
			io.db(),
			io.readOpt(),
			io.writeOpt(),
			io.registerObject()
		)
	),
	dict_
	(
		IOobject
		(
			searchableSurface::name() + "Dict",
			searchableSurface::instance(),
			searchableSurface::local(),
			searchableSurface::db(),
			searchableSurface::readOpt(),
			searchableSurface::writeOpt(),
			searchableSurface::registerObject()
		)
	)
{
	read();

	if (debug)
	{
		Info<< "Read distributedTriSurface from " << io.objectPath()
			<< ':' << endl;
		writeStats(Info);

		labelList nTris(Pstream::nProcs());
		nTris[Pstream::myProcNo()] = triSurface::size();
		Pstream::gatherList(nTris);
		Pstream::scatterList(nTris);

		Info<< endl<< "\tproc\ttris\tbb" << endl;
		forAll(nTris, procI)
		{
			Info<< '\t' << procI << '\t' << nTris[procI]
				 << '\t' << procBb_[procI] << endl;
		}
		Info<< endl;
	}
}


Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
(
	const IOobject& io,
	const dictionary& dict
)
:
	//triSurfaceMesh(io, dict),
	triSurfaceMesh
	(
		IOobject
		(
			io.name(),
			io.time().findInstance(io.local(), word::null),
			io.local(),
			io.db(),
			io.readOpt(),
			io.writeOpt(),
			io.registerObject()
		),
		dict
	),
	dict_
	(
		IOobject
		(
			searchableSurface::name() + "Dict",
			searchableSurface::instance(),
			searchableSurface::local(),
			searchableSurface::db(),
			searchableSurface::readOpt(),
			searchableSurface::writeOpt(),
			searchableSurface::registerObject()
		)
	)
{
	read();

	if (debug)
	{
		Info<< "Read distributedTriSurface from " << io.objectPath()
			<< " and dictionary:" << endl;
		writeStats(Info);

		labelList nTris(Pstream::nProcs());
		nTris[Pstream::myProcNo()] = triSurface::size();
		Pstream::gatherList(nTris);
		Pstream::scatterList(nTris);

		Info<< endl<< "\tproc\ttris\tbb" << endl;
		forAll(nTris, procI)
		{
			Info<< '\t' << procI << '\t' << nTris[procI]
				 << '\t' << procBb_[procI] << endl;
		}
		Info<< endl;
	}
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::distributedTriSurfaceMesh::~distributedTriSurfaceMesh()
{
	clearOut();
}


void Foam::distributedTriSurfaceMesh::clearOut()
{
	globalTris_.clear();
	triSurfaceMesh::clearOut();
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
{
	if (!globalTris_.valid())
	{
		globalTris_.reset(new globalIndex(triSurface::size()));
	}
	return globalTris_;
}


void Foam::distributedTriSurfaceMesh::findNearest
(
	const pointField& samples,
	const scalarField& nearestDistSqr,
	List<pointIndexHit>& info
) const
{
	const indexedOctree<treeDataTriSurface>& octree = tree();

	// Important:force synchronised construction of indexing
	const globalIndex& triIndexer = globalTris();


	// Initialise
	// ~~~~~~~~~~

	info.setSize(samples.size());
	forAll(info, i)
	{
		info[i].setMiss();
	}



	// Do any local queries
	// ~~~~~~~~~~~~~~~~~~~~

	label nLocal = 0;

	{
		// Work array - whether processor bb overlaps the bounding sphere.
		boolList procBbOverlaps(Pstream::nProcs());

		forAll(samples, i)
		{
			// Find the processor this sample+radius overlaps.
			label nProcs = calcOverlappingProcs
			(
				samples[i],
				nearestDistSqr[i],
				procBbOverlaps
			);

			// Overlaps local processor?
			if (procBbOverlaps[Pstream::myProcNo()])
			{
				info[i] = octree.findNearest(samples[i], nearestDistSqr[i]);
				if (info[i].hit())
				{
					info[i].setIndex(triIndexer.toGlobal(info[i].index()));
				}
				if (nProcs == 1)
				{
					// Fully local
					nLocal++;
				}
			}
		}
	}


	if
	(
		Pstream::parRun()
	 && (
			returnReduce(nLocal, sumOp<label>())
		  < returnReduce(samples.size(), sumOp<label>())
		)
	)
	{
		// Not all can be resolved locally. Build queries and map, send over
		// queries, do intersections, send back and merge.

		// Calculate queries and exchange map
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		pointField allCentres;
		scalarField allRadiusSqr;
		labelList allSegmentMap;
		autoPtr<mapDistribute> mapPtr
		(
			calcLocalQueries
			(
				samples,
				nearestDistSqr,

				allCentres,
				allRadiusSqr,
				allSegmentMap
			)
		);
		const mapDistribute& map = mapPtr();


		// swap samples to local processor
		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		map.distribute
		(
			//Pstream::scheduled,
			//map.schedule(),
			Pstream::nonBlocking,
			List<labelPair>(0),
			map.constructSize(),
			map.subMap(),           // what to send
			map.constructMap(),     // what to receive
			allCentres
		);
		map.distribute
		(
			//Pstream::scheduled,
			//map.schedule(),
			Pstream::nonBlocking,
			List<labelPair>(0),
			map.constructSize(),
			map.subMap(),           // what to send
			map.constructMap(),     // what to receive
			allRadiusSqr
		);


		// Do my tests
		// ~~~~~~~~~~~

		List<pointIndexHit> allInfo(allCentres.size());
		forAll(allInfo, i)
		{
			allInfo[i] = octree.findNearest
			(
				allCentres[i],
				allRadiusSqr[i]
			);
			if (allInfo[i].hit())
			{
				allInfo[i].setIndex(triIndexer.toGlobal(allInfo[i].index()));
			}
		}


		// Send back results
		// ~~~~~~~~~~~~~~~~~

		map.distribute
		(
			//Pstream::scheduled,
			//map.schedule            // note reverse schedule
			//(
			//    map.constructMap(),
			//    map.subMap()
			//),
			Pstream::nonBlocking,
			List<labelPair>(0),
			allSegmentMap.size(),
			map.constructMap(),     // what to send
			map.subMap(),           // what to receive
			allInfo
		);


		// Extract information
		// ~~~~~~~~~~~~~~~~~~~

		forAll(allInfo, i)
		{
			if (allInfo[i].hit())
			{
				label pointI = allSegmentMap[i];

				if (!info[pointI].hit())
				{
					// No intersection yet so take this one
					info[pointI] = allInfo[i];
				}
				else
				{
					// Nearest intersection
					if
					(
					    magSqr(allInfo[i].hitPoint()-samples[pointI])
					  < magSqr(info[pointI].hitPoint()-samples[pointI])
					)
					{
					    info[pointI] = allInfo[i];
					}
				}
			}
		}
	}
}


void Foam::distributedTriSurfaceMesh::findLine
(
	const pointField& start,
	const pointField& end,
	List<pointIndexHit>& info
) const
{
	findLine
	(
		true,   // nearestIntersection
		start,
		end,
		info
	);
}


void Foam::distributedTriSurfaceMesh::findLineAny
(
	const pointField& start,
	const pointField& end,
	List<pointIndexHit>& info
) const
{
	findLine
	(
		true,   // nearestIntersection
		start,
		end,
		info
	);
}


void Foam::distributedTriSurfaceMesh::findLineAll
(
	const pointField& start,
	const pointField& end,
	List<List<pointIndexHit> >& info
) const
{
	// Reuse fineLine. We could modify all of findLine to do multiple
	// intersections but this would mean a lot of data transferred so
	// for now we just find nearest intersection and retest from that
	// intersection onwards.

	// Work array.
	List<pointIndexHit> hitInfo(start.size());

	findLine
	(
		true,   // nearestIntersection
		start,
		end,
		hitInfo
	);

	// Tolerances:
	// To find all intersections we add a small vector to the last intersection
	// This is chosen such that
	// - it is significant (SMALL is smallest representative relative tolerance;
	//   we need something bigger since we're doing calculations)
	// - if the start-end vector is zero we still progress
	const vectorField dirVec(end-start);
	const scalarField magSqrDirVec(magSqr(dirVec));
	const vectorField smallVec
	(
		Foam::sqrt(SMALL)*dirVec
	  + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
	);

	// Copy to input and compact any hits
	labelList pointMap(start.size());
	pointField e0(start.size());
	pointField e1(start.size());
	label compactI = 0;

	info.setSize(hitInfo.size());
	forAll(hitInfo, pointI)
	{
		if (hitInfo[pointI].hit())
		{
			info[pointI].setSize(1);
			info[pointI][0] = hitInfo[pointI];

			point pt = hitInfo[pointI].hitPoint() + smallVec[pointI];

			if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
			{
				e0[compactI] = pt;
				e1[compactI] = end[pointI];
				pointMap[compactI] = pointI;
				compactI++;
			}
		}
		else
		{
			info[pointI].clear();
		}
	}

	e0.setSize(compactI);
	e1.setSize(compactI);
	pointMap.setSize(compactI);

	while (returnReduce(e0.size(), sumOp<label>()) > 0)
	{
		findLine
		(
			true,   // nearestIntersection
			e0,
			e1,
			hitInfo
		);


		// Extract
		label compactI = 0;
		forAll(hitInfo, i)
		{
			if (hitInfo[i].hit())
			{
				label pointI = pointMap[i];

				label sz = info[pointI].size();
				info[pointI].setSize(sz+1);
				info[pointI][sz] = hitInfo[i];

				point pt = hitInfo[i].hitPoint() + smallVec[pointI];

				if (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
				{
					e0[compactI] = pt;
					e1[compactI] = end[pointI];
					pointMap[compactI] = pointI;
					compactI++;
				}
			}
		}

		// Trim
		e0.setSize(compactI);
		e1.setSize(compactI);
		pointMap.setSize(compactI);
	}
}


void Foam::distributedTriSurfaceMesh::getRegion
(
	const List<pointIndexHit>& info,
	labelList& region
) const
{
	if (!Pstream::parRun())
	{
		region.setSize(info.size());
		forAll(info, i)
		{
			if (info[i].hit())
			{
				region[i] = triSurface::operator[](info[i].index()).region();
			}
			else
			{
				region[i] = -1;
			}
		}
		return;
	}


	// Get query data (= local index of triangle)
	// ~~~~~~~~~~~~~~

	labelList triangleIndex(info.size());
	autoPtr<mapDistribute> mapPtr
	(
		calcLocalQueries
		(
			info,
			triangleIndex
		)
	);
	const mapDistribute& map = mapPtr();


	// Do my tests
	// ~~~~~~~~~~~

	const triSurface& s = static_cast<const triSurface&>(*this);

	region.setSize(triangleIndex.size());

	forAll(triangleIndex, i)
	{
		label triI = triangleIndex[i];
		region[i] = s[triI].region();
	}


	// Send back results
	// ~~~~~~~~~~~~~~~~~

	map.distribute
	(
		//Pstream::scheduled,
		//map.schedule            // note reverse schedule
		//(
		//    map.constructMap(),
		//    map.subMap()
		//),
		Pstream::nonBlocking,
		List<labelPair>(0),
		info.size(),
		map.constructMap(),     // what to send
		map.subMap(),           // what to receive
		region
	);
}


void Foam::distributedTriSurfaceMesh::getNormal
(
	const List<pointIndexHit>& info,
	vectorField& normal
) const
{
	if (!Pstream::parRun())
	{
		triSurfaceMesh::getNormal(info, normal);
		return;
	}


	// Get query data (= local index of triangle)
	// ~~~~~~~~~~~~~~

	labelList triangleIndex(info.size());
	autoPtr<mapDistribute> mapPtr
	(
		calcLocalQueries
		(
			info,
			triangleIndex
		)
	);
	const mapDistribute& map = mapPtr();


	// Do my tests
	// ~~~~~~~~~~~

	const triSurface& s = static_cast<const triSurface&>(*this);

	normal.setSize(triangleIndex.size());

	forAll(triangleIndex, i)
	{
		label triI = triangleIndex[i];
		normal[i] = s[triI].normal(s.points());
		normal[i] /= mag(normal[i]) + VSMALL;
	}


	// Send back results
	// ~~~~~~~~~~~~~~~~~

	map.distribute
	(
		//Pstream::scheduled,
		//map.schedule            // note reverse schedule
		//(
		//    map.constructMap(),
		//    map.subMap()
		//),
		Pstream::nonBlocking,
		List<labelPair>(0),
		info.size(),
		map.constructMap(),     // what to send
		map.subMap(),           // what to receive
		normal
	);
}


void Foam::distributedTriSurfaceMesh::getField
(
	const List<pointIndexHit>& info,
	labelList& values
) const
{
	if (!Pstream::parRun())
	{
		triSurfaceMesh::getField(info, values);
		return;
	}

	if (foundObject<triSurfaceLabelField>("values"))
	{
		const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
		(
			"values"
		);


		// Get query data (= local index of triangle)
		// ~~~~~~~~~~~~~~

		labelList triangleIndex(info.size());
		autoPtr<mapDistribute> mapPtr
		(
			calcLocalQueries
			(
				info,
				triangleIndex
			)
		);
		const mapDistribute& map = mapPtr();


		// Do my tests
		// ~~~~~~~~~~~

		values.setSize(triangleIndex.size());

		forAll(triangleIndex, i)
		{
			label triI = triangleIndex[i];
			values[i] = fld[triI];
		}


		// Send back results
		// ~~~~~~~~~~~~~~~~~

		map.distribute
		(
			Pstream::nonBlocking,
			List<labelPair>(0),
			info.size(),
			map.constructMap(),     // what to send
			map.subMap(),           // what to receive
			values
		);
	}
}


void Foam::distributedTriSurfaceMesh::getVolumeType
(
	const pointField& points,
	List<volumeType>& volType
) const
{
	FatalErrorIn
	(
		"distributedTriSurfaceMesh::getVolumeType"
		"(const pointField&, List<volumeType>&) const"
	)   << "Volume type not supported for distributed surfaces."
		<< exit(FatalError);
}


// Subset the part of surface that is overlapping bb.
Foam::triSurface Foam::distributedTriSurfaceMesh::overlappingSurface
(
	const triSurface& s,
	const List<treeBoundBox>& bbs,

	labelList& subPointMap,
	labelList& subFaceMap
)
{
	// Determine what triangles to keep.
	boolList includedFace(s.size(), false);

	// Create a slightly larger bounding box.
	List<treeBoundBox> bbsX(bbs.size());
	const scalar eps = 1.0e-4;
	forAll(bbs, i)
	{
		const point mid = 0.5*(bbs[i].min() + bbs[i].max());
		const vector halfSpan = (1.0+eps)*(bbs[i].max() - mid);

		bbsX[i].min() = mid - halfSpan;
		bbsX[i].max() = mid + halfSpan;
	}

	forAll(s, triI)
	{
		const labelledTri& f = s[triI];
		const point& p0 = s.points()[f[0]];
		const point& p1 = s.points()[f[1]];
		const point& p2 = s.points()[f[2]];

		if (overlaps(bbsX, p0, p1, p2))
		{
			includedFace[triI] = true;
		}
	}

	return subsetMesh(s, includedFace, subPointMap, subFaceMap);
}


void Foam::distributedTriSurfaceMesh::distribute
(
	const List<treeBoundBox>& bbs,
	const bool keepNonLocal,
	autoPtr<mapDistribute>& faceMap,
	autoPtr<mapDistribute>& pointMap
)
{
	// Get bbs of all domains
	// ~~~~~~~~~~~~~~~~~~~~~~

	{
		List<List<treeBoundBox> > newProcBb(Pstream::nProcs());

		switch(distType_)
		{
			case FOLLOW:
				newProcBb[Pstream::myProcNo()].setSize(bbs.size());
				forAll(bbs, i)
				{
					newProcBb[Pstream::myProcNo()][i] = bbs[i];
				}
				Pstream::gatherList(newProcBb);
				Pstream::scatterList(newProcBb);
			break;

			case INDEPENDENT:
				newProcBb = independentlyDistributedBbs(*this);
			break;

			case FROZEN:
				return;
			break;

			default:
				FatalErrorIn("distributedTriSurfaceMesh::distribute(..)")
					<< "Unsupported distribution type." << exit(FatalError);
			break;
		}

		//if (debug)
		//{
		//    Info<< "old bb:" << procBb_ << endl << endl;
		//    Info<< "new bb:" << newProcBb << endl << endl;
		//    Info<< "Same:" << (newProcBb == procBb_) << endl;
		//}

		if (newProcBb == procBb_)
		{
			return;
		}
		else
		{
			procBb_.transfer(newProcBb);
			dict_.set("bounds", procBb_[Pstream::myProcNo()]);
		}
	}


	// Debug information
	if (debug)
	{
		labelList nTris(Pstream::nProcs());
		nTris[Pstream::myProcNo()] = triSurface::size();
		Pstream::gatherList(nTris);
		Pstream::scatterList(nTris);

		Info<< "distributedTriSurfaceMesh::distribute : before distribution:"
			<< endl
			<< "\tproc\ttris" << endl;

		forAll(nTris, procI)
		{
			Info<< '\t' << procI << '\t' << nTris[procI] << endl;
		}
		Info<< endl;
	}


	// Use procBbs to determine which faces go where
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	labelListList faceSendMap(Pstream::nProcs());
	labelListList pointSendMap(Pstream::nProcs());

	forAll(procBb_, procI)
	{
		overlappingSurface
		(
			*this,
			procBb_[procI],
			pointSendMap[procI],
			faceSendMap[procI]
		);

		if (debug)
		{
			//Pout<< "Overlapping with proc " << procI
			//    << " faces:" << faceSendMap[procI].size()
			//    << " points:" << pointSendMap[procI].size() << endl << endl;
		}
	}

	if (keepNonLocal)
	{
		// Include in faceSendMap/pointSendMap the triangles that are
		// not mapped to any processor so they stay local.

		const triSurface& s = static_cast<const triSurface&>(*this);

		boolList includedFace(s.size(), true);

		forAll(faceSendMap, procI)
		{
			if (procI != Pstream::myProcNo())
			{
				forAll(faceSendMap[procI], i)
				{
					includedFace[faceSendMap[procI][i]] = false;
				}
			}
		}

		// Combine includedFace (all triangles that are not on any neighbour)
		// with overlap.

		forAll(faceSendMap[Pstream::myProcNo()], i)
		{
			includedFace[faceSendMap[Pstream::myProcNo()][i]] = true;
		}

		subsetMesh
		(
			s,
			includedFace,
			pointSendMap[Pstream::myProcNo()],
			faceSendMap[Pstream::myProcNo()]
		);
	}


	// Send over how many faces/points I need to receive
	// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	labelListList faceSendSizes(Pstream::nProcs());
	faceSendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
	forAll(faceSendMap, procI)
	{
		faceSendSizes[Pstream::myProcNo()][procI] = faceSendMap[procI].size();
	}
	Pstream::gatherList(faceSendSizes);
	Pstream::scatterList(faceSendSizes);



	// Exchange surfaces
	// ~~~~~~~~~~~~~~~~~

	// Storage for resulting surface
	List<labelledTri> allTris;
	pointField allPoints;

	labelListList faceConstructMap(Pstream::nProcs());
	labelListList pointConstructMap(Pstream::nProcs());


	// My own surface first
	// ~~~~~~~~~~~~~~~~~~~~

	{
		labelList pointMap;
		triSurface subSurface
		(
			subsetMesh
			(
				*this,
				faceSendMap[Pstream::myProcNo()],
				pointMap
			)
		);

		allTris = subSurface;
		allPoints = subSurface.points();

		faceConstructMap[Pstream::myProcNo()] = identity
		(
			faceSendMap[Pstream::myProcNo()].size()
		);
		pointConstructMap[Pstream::myProcNo()] = identity
		(
			pointSendMap[Pstream::myProcNo()].size()
		);
	}



	// Send all
	// ~~~~~~~~

	forAll(faceSendSizes, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			if (faceSendSizes[Pstream::myProcNo()][procI] > 0)
			{
				OPstream str(Pstream::blocking, procI);

				labelList pointMap;
				triSurface subSurface
				(
					subsetMesh
					(
					    *this,
					    faceSendMap[procI],
					    pointMap
					)
				);

				//if (debug)
				//{
				//    Pout<< "Sending to " << procI
				//        << " faces:" << faceSendMap[procI].size()
				//        << " points:" << subSurface.points().size() << endl
				//        << endl;
				//}

				str << subSurface;
		   }
		}
	}


	// Receive and merge all
	// ~~~~~~~~~~~~~~~~~~~~~

	forAll(faceSendSizes, procI)
	{
		if (procI != Pstream::myProcNo())
		{
			if (faceSendSizes[procI][Pstream::myProcNo()] > 0)
			{
				IPstream str(Pstream::blocking, procI);

				// Receive
				triSurface subSurface(str);

				//if (debug)
				//{
				//    Pout<< "Received from " << procI
				//        << " faces:" << subSurface.size()
				//        << " points:" << subSurface.points().size() << endl
				//        << endl;
				//}

				// Merge into allSurf
				merge
				(
					mergeDist_,
					subSurface,
					subSurface.points(),

					allTris,
					allPoints,
					faceConstructMap[procI],
					pointConstructMap[procI]
				);

				//if (debug)
				//{
				//    Pout<< "Current merged surface : faces:" << allTris.size()
				//        << " points:" << allPoints.size() << endl << endl;
				//}
		   }
		}
	}


	faceMap.reset
	(
		new mapDistribute
		(
			allTris.size(),
			faceSendMap,
			faceConstructMap,
			true
		)
	);
	pointMap.reset
	(
		new mapDistribute
		(
			allPoints.size(),
			pointSendMap,
			pointConstructMap,
			true
		)
	);

	// Construct triSurface. Reuse storage.
	triSurface::operator=(triSurface(allTris, patches(), allPoints, true));

	clearOut();

	// Regions stays same
	// Volume type stays same.

	distributeFields<label>(faceMap());
	distributeFields<scalar>(faceMap());
	distributeFields<vector>(faceMap());
	distributeFields<sphericalTensor>(faceMap());
	distributeFields<symmTensor>(faceMap());
	distributeFields<tensor>(faceMap());

	if (debug)
	{
		labelList nTris(Pstream::nProcs());
		nTris[Pstream::myProcNo()] = triSurface::size();
		Pstream::gatherList(nTris);
		Pstream::scatterList(nTris);

		Info<< "distributedTriSurfaceMesh::distribute : after distribution:"
			<< endl
			<< "\tproc\ttris" << endl;

		forAll(nTris, procI)
		{
			Info<< '\t' << procI << '\t' << nTris[procI] << endl;
		}
		Info<< endl;
	}
}


//- Write using given format, version and compression
bool Foam::distributedTriSurfaceMesh::writeObject
(
	IOstream::streamFormat fmt,
	IOstream::versionNumber ver,
	IOstream::compressionType cmp
) const
{
	// Make sure dictionary goes to same directory as surface
	const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();

	// Dictionary needs to be written in ascii - binary output not supported.
	return
		triSurfaceMesh::writeObject(fmt, ver, cmp)
	 && dict_.writeObject(IOstream::ASCII, ver, cmp);
}


void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
{
	boundBox bb;
	label nPoints;
	calcBounds(bb, nPoints);
	reduce(bb.min(), minOp<point>());
	reduce(bb.max(), maxOp<point>());

	os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
		<< endl
		<< "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
		<< "Bounding Box : " << bb << endl;
}


// ************************************************************************* //
